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1.  Overview  and  Project  Goals 


Our  research  investigated  the  connections  between  scale  space  and  the  linear  and 
nonlinear  diffusion  of  images  using  a  combination  of  analytic  methods  and  prototype 
Matlab  and  Mathematica  programs.  We  investigated  concepts  of  generalized  entropy  and 
quantum  entropy  and  their  measures,  using  mathematical  analysis.  The  generalized 
entropies  included  in  particular  Renyi  entropy  [renyi]  with  parameter  q.  We  also 
examined  the  properties  of  generalized  entropy  of  images  subject  to  linear  and  nonlinear 
diffusion  processes.  We  performed  analysis  of  quantum  and  semi-classical  entropies  of 
model  physical  systems. 

We  investigated  the  feasibility  of  applying  forms  of  generalized  quantum  search  to 
scheduling  and  logistics  problems.  As  part  of  this  effort,  we  performed  simulations  of 
adaptive  quantum  search.  Based  upon  this  portion  of  the  research,  Ben  J.  Jones  received 
the  Ryan  Sayers  memorial  award  for  research  as  a  senior  student  majoring  in  both 
computer  science  and  engineering  physics. 

As  a  further  goal,  we  investigated  the  capabilities  and  constraints  of  quantum  lattice  gas 
algorithms  (QLGAs),  including  the  properties  of  their  (unitary)  collision  operators. 
QLGAs  are  known  to  be  able  to  simulate  a  variety  of  partial  differential  equations 
[yepez],  including  the  Navier-Stokes,  Boltzmann,  and  Dirac  equations  [iwoprd].  In  our 
research,  we  investigated  the  feasibility  of  the  solution  of  the  celebrated  Maxwell 
equations  of  electromagnetism  by  means  of  QLGAs  [coffey].  We  also  formulated, 
implemented,  and  simulated  a  QLGA  for  the  telegraph  equation  with  one  spatial 
dimension  [cc2009].  The  latter  equation  combines  facets  of  both  the  wave  and  diffusion 
equations.  Our  research  indicated  that  the  simulation  of  the  ID  linear  telegraph  equation 
with  an  iterate  constraint  has  application  to  signal  denoising. 
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2.  Approach 


Consideration  of  classical  and  quantum  entropy  and  their  application  motivated  much  of 
our  initial  effort.  Furthermore,  entropy  concepts  occurred  throughout  much  of  the  later 
research.  Most  explicitly  they  appear  in  the  consideration  of  image  entropy,  but  also  later 
for  instance  in  the  Schmidt  strength  from  quantum  logic  gate  decomposition.  This  form 
of  entropy  gives  a  measure  of  the  nonlocal  content  of  an  entangling  logic  gate. 

Entropy  and  its  generalizations  are  among  the  most  important  measures  of  information 
content  and  complexity  of  signals  and  images.  Another  important  concept  that  we  made 
use  of  was  that  of  scale  spaces.  In  scale-space  theory  one  embeds  an  image  into  a 
continuous  family  of  gradually  smoother  versions  of  it.  The  time  t  acts  as  a  parameter  for 
this,  with  the  original  image  corresponding  to  t  =  0.  Increasing  the  scale  should  simplify 
the  image  without  creating  spurious  structures.  For  instance,  in  viewing  a  facial  image  at 
coarser  scales,  it  would  be  undesirable  to  have  artificial  features  appearing.  A  scale- 
space  introduces  a  hierarchy  of  image  features,  and  can  provide  an  important  process  in 
going  from  a  pixel-level  description  to  a  semantical  image  description. 

Partial  differential  equations  (PDEs)  are  the  suitable  framework  for  scale-spaces,  and  the 
oldest,  simplest,  and  probably  most  studied  version  of  scale  space  corresponds  to  a  linear 
diffusion  process.  The  fundamental  solution  (Greens  function)  for  a  linear  heat  or 
diffusion  equation  is  a  Gaussian  function  with  standard  deviation  proportional  to  the 
square  root  of  the  time.  The  solution  of  the  linear  diffusion  equation  can  be  given  as  the 
convolution  of  the  initial  data  (image)  with  this  Gaussian  function,  and  this  gives  linear 
scale  space.  For  the  extension  to  nonlinear  scale  space,  the  PDE  involved  is  a  nonlinear 
diffusion  equation  with  a  decreasing  diffusion  coefficient  D.  When  D  depends  upon  the 
magnitude  of  the  gradient  of  the  pixel  intensities,  we  generally  obtain  selective 
smoothing.  Eocations  where  the  gradient  is  large  have  strong  probability  of  being  an 
edge,  and  D  is  reduced. 

A  form  of  the  diffusion  (or  heat)  equation  is  dVdi  =  V  •  [D(I)  •  VI],  where  I  is  the  image 
intensity.  For  linear  and  isotropic  diffusion  in  two  dimensions,  the  diffusion  equation  is 
simply  It  =  D(Ixx+Iyy).  The  basic  idea  of  anisotropic  diffusion  is  to  diffuse  intensities 
along  edges  of  objects  that  appear  within  an  image,  while  not  diffusing  (or  even 
enhancing  the  contrast)  along  directions  that  are  perpendicular  to  edges  [perona,  price]. 
For  nonlinear  diffusion  of  images  to  generate  intraregion  smoothing,  D  =  D(|VI|),  and 
regions  of  high  intensity  contrast  undergo  less  diffusion. 

The  diffusion  boundary  value  problem  is  completed  with  the  specification  of  initial  and 
boundary  conditions.  The  initial  condition  is  simply  l(x,y,t=0)  =  f(x,y),  where  f  is  the 
given  image.  Homogeneous  Neumann  or  periodic  boundary  conditions  are  appropriate  as 
they  lead  to  conservation  of  the  average  grey  values  of  the  whole  image.  This  fact  results 
from  the  divergence  form  of  the  diffusion  equation. 
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In  the  mathematical  sense,  the  partial  differential  equations  we  simulated  with  QLGAs 
are  evolution  equations  of  the  form 


Ut  =  f(u,Ux,Uxx),  (1) 

where  the  subscripts  denote  partial  differentiation,  for  example,  Ut  =  dvJdi.  Therefore,  we 
solved  equations  that  are  first  order  in  time  and  second  order  in  space.  The  PDE  is 
supplemented  by  an  initial  condition  u(x,t  =  0)  =  u(x,0)  and  boundary  conditions  at  the 
endpoints  of  the  x-interval  to  completely  specify  the  problem.  We  used  most  often 
periodic  boundary  conditions,  so  that  for  a  one-dimensional  problem  on  an  interval  of 
length  L,  u(0,t)  =  u(L,t).  For  image  processing  applications  in  particular,  the  dependent 
variable  u  corresponds  to  pixel  intensities.  For  heat  transfer  problems,  u  represents  the 
temperature,  while  for  fluid  flow  problems  u  corresponds  to  either  the  mass  density  or 
flow  velocity. 

A  quantum  lattice  gas  algorithm  is  a  quantum  version  of  a  classical  lattice  gas,  which  in 
turn  is  an  extension  of  classical  cellular  automata  [yepez,  doolen].  In  place  of  the  binary 
lattice  variables  of  a  classical  lattice  gas,  the  quantum  version  has  a  local  Hilbert  space 
describing  the  quantum  bit  (qubits).  In  the  classical  case,  in  order  to  recover  the  proper 
macroscopic  dynamics  (and  thermodynamics),  it  is  important  to  ensure  that  the 
microscopic  dynamics  preserves  conservation  laws.  Similarly,  in  the  quantum  case,  the 
number  densities  of  qubits  must  be  preserved  so  that  the  interaction  operator,  called  the 
collision  operator,  must  be  unitary. 

The  operation  of  a  type-II  quantum  processor  includes  the  sequential  repetition  of  four 
main  steps  [berman,  yepez,  vahala,  love].  First,  initialization  creates  the  quantum- 
mechanical  initial  state  that  corresponds  to  the  initial  probability  distribution  for  a  partial 
differential  equation  to  be  solved.  Secondly,  a  unitary  transformation  (collision  operator) 
is  applied  in  parallel  to  all  the  local  Hilbert  spaces  in  the  lattice.  Next,  in  the 
measurement  step,  the  quantum  states  of  all  the  nodes  are  read  out.  Fastly,  these  results 
are  streamed  to  neighboring  sites  to  reinitialize  the  quantum  processor  in  the  state  which 
corresponds  to  the  new  probability  distribution. 

After  the  initialization  step,  a  quantum  lattice  gas  algorithm  performs  iterations  of  a 
collision  operator  and  a  streaming  operator.  The  latter  operator  shifts  the  state  of  a  qubit 
from  a  given  lattice  site  to  its  nearest  neighbors  in  the  lattice. 
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3.  Details  of  Technical  Accomplishments 

3.1  Entropy  investigation 

3.1.1  Diffusion  processing 

Diffusion  processing  has  proven  very  useful  for  praetieal  image  enhaneement,  wherein 
the  visual  quality  of  an  image  is  improved.  We  have  investigated  methods  of  carrying 
out  such  processing  in  a  combined  elassical-quantum  eomputing  environment. 

We  performed  investigations  speeifically  examining  the  ehange  in  generalized  entropies, 
including  Renyi  entropy  [renyi].  We  are  motivated  to  examine  the  Renyi  entropy 
because  it  seems  to  comport  well  with  the  ideas  of  seale  spaees  diseussed  above.  Like 
Tsallis  entropy,  Renyi  entropy  eontains  a  parameter  q  1;  Sq(p)  =  ln[Sj=i’^  Pj‘*]/(l-q), 
where  pj  are  normalized  pixel  values.  It  also  reduees  to  the  Shannon- Wiener  entropy 
when  q  ^  1:  Ssw(p)  =  Pj  In  pj.  Here  we  have  sums  over  pixel  values  (probabilities) 
sinee  we  are  in  the  diserete  (digital)  ease  rather  than  integrals  suitable  to  the  eontinuous 
ease.  We  see  that  the  Renyi  entropy  gives  a  sort  of  extrapolation  of  the  Shannon- Wiener 
entropy  [karolz].  The  parameter  q  serves  to  introduee  a  type  of  weighting  of  pixel  values 
into  the  entropy  and  may  help  to  provide  a  hierarehical  ordering  of  the  image  information. 
Therefore,  it  is  directly  connected  to  seale  spaees. 

Figure  1  shows  an  example  of  the  deerease  of  various  entropies  upon  sueeessive 
iterations  of  nonlinear  diffusion  with  an  exponential  diffusivity  funetion.  The  image 
entropies  are  determined  from  the  histogram  of  the  distribution  of  pixel  intensities. 
Matlab’s  entropy  just  uses  a  different  base  (2)  of  logarithm  than  Shannon.  The  ehange  in 
image  entropy  during  diffusion  proeessing  is  a  refleetion  of  the  loss  of  information  in 
going  from  finer  to  eoarser  scales  of  resolution. 

We  may  list  the  following  as  desirable  properties  of  entropy  as  applied  to  image 
processing  tasks.  These  points  could  serve  as  guidelines  for  future  researeh  (cf  [starek]), 
or  as  a  start  on  axioms  for  developing  measures  of  information  content. 

•  A  single-value  image  has  zero  information  eontent. 

•  The  amount  of  information  in  an  image  is  independent  of  the 
baekground. 

•  The  amount  of  information  is  dependent  on  the  noise. 

•  The  entropy  should  work  in  the  same  way  for  a  pixel  that  has  value  V  +  5  and 
for  a  pixel  with  value  V  -  5,  with  V  being  the  background  value. 

•  The  amount  of  information  is  dependent  on  the  spatial  eorrelation  in  the 
image.  An  image  with  large  homogeneous  features  (above  the  noise) 
eontains  less  information. 
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Figure  1.  Decrease  of  generalized  image  entropy  with  iteration  nnmber  dnring  nonlinear  diffnsion 

processing. 
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3.1.2  Entropy  calculations  for  model  physical  systems 

Lately  information  theoretie  eoneepts  have  been  playing  a  larger  role  in  quantum 
meehanies  and  quantum  eomputing.  Information  eoneepts  have  been  employed  in  both 
fundamental  diseussions  and  in  praetieal  applieations  ineluding  synthesis  and  analysis  of 
eleetron  densities  in  position  and  momentum  spaee  [gadre].  Indeed,  the  sum  of  quantum 
position  and  momentum  entropies  has  been  advoeated  as  a  measure  of  wavefunetion 
quality  [gadre]. 

Quantum  entropy  provides  a  quantitative  deseription  about  the  uneertainty  or  laek  of 
knowledge  of  an  observable.  For  instanee,  the  entropy  is  one  measure  of  the 
deloealization  of  a  wavepaeket.  The  Shannon  entropy  provides  an  unambiguous  measure 
whieh  is  eomplementary  to  the  information  eontent  of  a  system.  With  the  reeent  researeh 
into  quantum  eomputing  [kitaev,  nielsen,  qereviews]  there  has  been  added  attention  on 
the  fundamental  physieal  limits  of  eomputation  [eofpla]  and  here  also  the  quantum 
entropy  plays  a  role. 

We  reeall  that  Hirsehman  [hirsehman]  antieipated  a  strengthened  quantum  uneertainty 
prineiple.  Later,  Deutseh  [deutseh]  and  others  [bbm]  showed  that  the  Heisenberg 
inequality  does  not  properly  express  the  quantum  uneertainty  prineiple  and  is  generally 
too  weak.  They  introdueed  entropy  measures  sueh  as  we  use  for  noneommuting 
observables. 

In  [eoffeyejp]  we  established  the  semielassieal  position  and  momentum  information 
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entropies  for  a  family  of  systems  with  rational  potential  energies  and  for  the  seeh 
potential  energy.  The  latter  is  an  important  instanee  sinee  it  applies  to  a  potential  with 
nonpolynomial  form.  It  also  ineludes  the  ease  of  an  attraetive  delta-funetion  potential  for 
eertain  limit  values  of  the  potential  parameters  of  strength  and  width.  The  resulting 
semielassieal  entropy  relations  have  high  utility.  This  is  beeause  otherwise  numerieal 
eomputation  may  be  required,  or  when  elosed  form  results  for  quantum  systems  are 
available,  the  multiple  sum  and  produet  expressions  do  not  readily  yield  physieal 
information. 

A  further  motivation  of  our  investigation  was  a  reeent  presentation  of  the  ground  state 
position  entropy  of  the  Poesehl-Teller  potential  [atre].  The  hyperbolie  form  of  this 
potential  is  none  other  than  the  seeh  funetional  dependenee  that  we  employ.  The  study 
[atre]  gave  some  numerieal  results  for  exeited  states,  although  we  note  that  the  exaet 
general  exeited  state  solution  is  expressible  in  terms  of  produets  of  operators  Op  =  d/dx  - 
p  tanh  ax  [lamb].  Sinee  asymptotie  relations  are  now  known  eonneeting  quantum 
entropies  to  elassieal  eounterparts  for  both  position  [ruiz]  and  momentum 
[eofmomentum],  we  are  able  to  derive  semielassieal  expressions  applieable  for  high 
exeited  states. 
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The  entropic  uncertainty  relation  gives  a  lower  bound  to  the  sum  of  position  and 
momentum  entropies  [bbm], 

Sq^’^4Sq^p^>D(1 +ln7r),  (2) 

for  a  D-dimensional  system.  This  inequality  stresses  the  reciprocity  of  position  and 
momentum  spaces.  For  if  the  wave  function  is  concentrated  in  coordinate  space,  it  will 
necessarily  be  more  diffuse  in  momentum  space,  and  vice  versa.  Not  surprisingly,  the 
bound  in  Equation  2  is  attained  by  Gaussian  wave  functions  [bbm].  Relation  2  extends  to 
other  pairs  of  noncommuting  observables  A  and  B,  S(A)  +  §(B)  > 

Sab,  where  Sab  is  a 

positive  constant.  From  Equation  2  follows  the  Heisenberg  uncertainty  relation,  showing 
that  this  inequality  is  stronger. 

The  entropy  sum  appears  in  many  other  lower  and  upper  bounds.  Another  lower  bound 
including  measurement  device  resolutions  Ax  of  position  and  Apx  of  momentum  is 
[bbirula] 

>  1  -  In  2  -In  Ax  Apx.  (3) 

An  upper  bound  to  the  entropy  sum  can  be  prescribed  in  terms  of  the  second  moments  in 
position  and  momentum  space  [gadre].  We  note  that  the  entropic  uncertainty  relation  has 
been  extended  to  nonzero  temperatures  T  [abesuzuki]. 

In  the  case  of  the  sech  potential  we  recall  that  the  classical  period  of  motion  depends 
upon  total  energy  E  as  T(E)  |E|'  .  We  then  determine  how  the  classical  position 
entropy  varies  logarithmically  with  the  energy.  By  invoking  an  asymptotic  relation 
[ruiz],  we  then  know  how  the  quantum  entropy  of  position  varies  for  the  nth  energy 
eigenstate.  Since  we  also  know  the  values  of  the  energy  levels  En,  we  determine 
explicitly  the  dependence  of  the  semiclassical  position  entropy  upon  principle  quantum 
number  n.  We  find  the  momentum  entropy  and  therefore  for  highly  excited 
states  [cofmomentum],  to  contain  In  |E|'^^^  dependence,  which  is  not  surprising  given  that 
the  Hamiltonian  is  quadratic  in  momentum  p.  In  the  paper  [coffeycjp]  all  these  relations 
are  made  quantitative.  The  calculations  require  fairly  advanced  integration  techniques. 
In  particular,  we  extended  several  known  results  employing  the  Gauss  hypergeometric 
function. 

Semiclassical  results  are  also  possible  and  discussed  for  a  family  of  rational  potentials  in 
[coffeycjp].  The  treatment  of  this  set  of  potentials  also  relies  heavily  upon  the  analytic 
properties  of  the  hypergeometric  function.  This  family  contains  a  parameter  q,  V(x)  =  - 
Vo/[l+(|x|/a)‘*],  where  Vo>  0,  such  that  the  limit  q  ^  oo  yields  the  finite  square  well  as 
shown  in  Figure  2.  For  q  =  2  the  rational  potential  is  very  similar  to  the  sech^  potential 
about  the  common  minimum  at  x  =  0. 
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x/a 

Figure  2.  The  sequence  of  potentials  V(x)A^o=  -l/[l+(|x|/a)‘']  obtained  for  varying  values  of  q  is 

plotted  versus  x/a. 

3,2  New  quantum  algorithm  for  the  modal  value 

In  this  section  we  present  a  quantum  algorithm  for  finding  the  most  often  occurring  (or 
modal)  value  of  a  data  set.  We  thereby  supplement  other  algorithms  that  can  determine 
the  mean  value  or  similar  quantities.  Our  algorithm  requires  the  combined  use  of 
quantum  counting  and  extended  quantum  search. 

The  mode  is  the  most  often  occurring  value  in  a  data  set  and  is  an  important  statistic.  For 
the  sake  of  definiteness  in  the  description,  we  assume  a  data  list  of  N  elements,  each  entry 
being  an  integer  in  the  range  [l,d].  We  propose  an  algorithm  that  uses  a  combination  of 
quantum  counting  [brassard]  and  quantum  search  [boyer,  groverl,  grover2,  groverS, 
nielsen],  and  makes  use  of  the  result  that  quantum  search  may  be  applied  to  find  more 
than  one  target  item  in  an  unsorted  list  and  that  the  number  of  target  items  need  not  be 
known  beforehand.  Our  method  gives  an  operational  complexity  of  O(dVN),  as  measured 
in  the  number  of  oracle  calls.  The  modal  value  need  not  be  unique  for  our  algorithm  to 
succeed. 
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Our  algorithm  for  determining  the  mode  is  related  to  other  algorithms  that  also  apply 
quantum  seareh,  give  a  quadratie  speed  up  over  the  elassieal  situation,  and  deliver  other 
useful  statisties.  These  inelude  an  algorithm  for  finding  the  minimum  or  maximum  value 
in  a  data  set  [durr].  In  addition,  the  mean  value  may  be  determined  by  [grover2],  and 
then  the  varianee  and  higher  moments  of  a  data  set  may  be  determined.  This  is  beeause 
the  moments  may  be  appropriately  written  as  averages.  In  partieular,  the  varianee  is  the 
mean  value  of  the  quantity  (x-<x>)^,  where  <x>  is  the  mean  value.  The  median  value 
may  also  be  estimated,  with  an  algorithm  giving  a  nearly  quadratie  speed  up  over 
elassieal  algorithms  in  the  worst  ease  [nayakwu]. 

Beeause  of  the  utility  to  eompute  averages,  quantum  seareh  also  has  applieation  to 
integration.  On  the  other  hand,  we  are  interested  to  have  quantum  algorithms  to 
determine  a  large  variety  of  statisties  of  given  data.  Sueh  statisties  eould  be  used  for 
instanee  in  analyzing  digital  images  by  applying  them  to  the  pixel  intensities.  Operations 
sueh  as  thresholding  and  region  segmentation  or  dilation  and  erosion  eould  be  performed 
based  upon  the  statistieal  values. 

We  reeall  that  quantum  eounting  relies  on  phase  estimation  and  the  Fourier  transform. 
Quantum  eounting  exploits  the  periodieity  with  the  number  of  Grover  iterates  of  the 
probability  amplitude  of  the  target  state(s).  In  turn,  quantum  phase  estimation  makes  use 
of  eontrolled-U^^  operations  and  the  inverse  Fourier  transform  to  give  the  best  n-bit 
estimate  of  the  phase  (j)  of  the  eigenvalue  e“'’  of  the  unitary  operator  U. 

The  quantum  seareh  algorithm  has  been  extended  in  a  number  of  ways,  ineluding  with 
different  iteration  operators  and  different  seleetive  phase  shifts  (e.g.,  [biham,  galindo, 
grover2]).  These  are  implementation  speeifie  matters  that  are  not  the  foeus  of  this 
diseussion. 


3.2.1  The  algorithm 

For  simplieity  of  deseription,  we  assume  that  no  value,  ineluding  speeifieally  the  modal 
value,  oeeurs  more  than  N/2  times  in  the  data  set.  This  is  not  a  limitation,  sinee  if  this 
number  M  >  N/2,  the  number  of  items  in  the  data  set  may  be  doubled  with  N  non-solution 
elements,  and  a  modified  oraele  for  quantum  seareh  may  be  suitably  eonstrueted  [nielsen]. 

We  assume  that  the  data  values  are  bounded  by  the  range  [l,d].  We  eonsider  two 
quantum  registers  |D)  and  |C),  the  data  and  eount  registers,  of  size  0(log2  N)  and  0(log2  d) 
qubits  respeetively.  The  algorithm  has  two  main  steps. 

First,  we  initialize  the  quantum  system  as  |D)|0).  Then  for  eaeh  data  value  in  the  first 
register  we  use  quantum  eounting  [brassard]  to  determine  if  that  value  oeeurs  and  if  so 
how  many  times,  sueh  eount  being  kept  in  the  seeond  register.  We  then  apply  the 
quantum  algorithm  for  determining  the  maximum  value  [durr]  that  uses  the  generalized 
quantum  seareh  algorithm  [boyer]  to  the  |C)  register.  This  returns  the  number  of  times  M 
the  modal  value  oeeurs. 
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Secondly,  we  reinitialize  the  system  as  |D)|0).  We  then  execute  an  extended  quantum 
search  to  find  the  data  element(s)  corresponding  to  the  previously  determined  number  of 
times  that  the  modal  value  occurs.  In  this  case,  the  oracle  function  indicates  which  target 
item(s)  occurs  M  times. 

This  algorithm  delivers  both  M  and  the  corresponding  modal  value(s)  of  the  data  set.  By 
the  use  of  extended  quantum  search,  the  modal  value  need  not  be  unique. 


3.2.2  Discussion  and  Summary 

The  probability  of  success  of  our  algorithm  can  be  boosted  to  be  arbitrarily  close  to  1 
throughout.  Suppose  we  desire  a  small  probability  of  failure  1  »  8  >  0.  To  guarantee 
that  all  of  our  counting  operations  succeed  with  probability  of  at  least  1-5,  we  require 
complexity  0(dA/Nlog(d/5)).  This  gives  our  algorithm  total  probability  >  (l/2)(l-5)  of 
success  [furrow].  From  here  we  appeal  to  amplitude  amplification  [bcwz]  for  maximum 
finding,  which  leads  to  an  overall  running  time  of  0(dA/Nlog(d/8)A/log(l/8},  where  8  is 
our  total  probability  of  failure. 

Our  work  serves  to  enlarge  the  collection  of  quantum  algorithms  for  determining 
statistics  of  a  data  set  to  include  the  modal  value.  A  combination  of  both  quantum 
counting  and  quantum  search  is  required  in  order  to  avoid  reducing  to  a  classical  0(N) 
complexity.  Our  description  here  is  not  intended  to  be  exhaustive  and  it  seems  likely  that 
several  variations  on  the  ideas  presented  are  possible. 


3,3  Quantum  logic  gate  decomposition 

Here  we  consider  two-qubit  operators  and  provide  a  correspondence  between  their 
Schmidt  number  and  controlled-NOT  (CNOT)  complexity,  where  the  CNOT  complexity 
is  up  to  local  unitary  operations.  The  results  are  obtained  by  complementary  means,  and 
a  number  of  examples  are  given.  For  full  details,  see  [coffeydeiotte]. 

An  operator  Q  acting  on  systems  A  and  B  may  be  written  as  the  operator-Schmidt 
decomposition 

Q  =  Ej  Sj  Aj  0  Bj,  (4) 

where  Sj  >  0  and  Aj  and  Bj  are  orthonormal  operator  bases  for  A  and  B,  respectively.  This 
form  may  be  proved  constructively  by  using  the  singular  value  decomposition.  Given  the 
representation  (4),  the  Schmidt  number  Sch(Q)  of  the  operator  Q  is  defined  as  the  number 
of  nonzero  coefficients  Sj. 

In  this  discussion  we  are  concerned  with  two-qubit  operators,  and  the  relation  of  their 
Schmidt  number  to  their  controlled-NOT  complexity.  It  is  known  that  two-qubit 
unitaries  are  equivalent  to  either  0,  1,  2,  or  3  CNOT  gates,  where  the  equivalence  is  up  to 
single-qubit  rotations  (e.g.,  [makhlin02]).  On  the  other  hand,  two-qubit  unitary  operators 
may  have  Schmidt  numbers  1,  2,  and  4,  but  not  3  [nielsenetal03].  (A  similar  result  for 
states  has  been  obtained  in  [dur02].) 
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We  recall  that  the  Schmidt  number  is  an  indicator  of  entanglement,  but  not  a  measure  of 
entanglement.  For  instance,  let  us  compare  the  Bell  state 

|(^)  =  (lW2)(|00)+|ll)), 

with  the  state 

|M/)  =  V(1-8')|00)  +  8|11), 

where  we  take  8  «  1 .  Both  of  these  states  have  Schmidt  number  two,  but  the  Bell  state 
is  much  more  entangled,  in  fact  maximally  entangled.  This  simple  example  reflects  a 
more  general  situation  where  a  single  term  dominates  in  the  Schmidt  decomposition.  The 
state  is  entangled,  but  it  may  be  weakly  so. 

In  [coffeydeiotte]  we  provide  a  classification  of  CNOT  complexity  viz  a  viz  the  operator- 
Schmidt  decomposition,  and  prove  it  in  alternative  ways.  We  then  provide  several 
examples,  and  discuss  related  topics.  Among  these,  we  present  relations  between 
different  approaches  for  finding  the  parameters  of  the  canonical  decomposition  of  two- 
qubit  operators. 

One  may  expect  that  as  the  Schmidt  number  increases,  so  too  does  the  CNOT  complexity. 
The  main  result  is  the  following  [coffeydeiotte].  Proposition:  A  two-qubit  operator  U 
with  CNOT  complexity  2  or  3  has  Sch(U)  =  4,  with  CNOT  complexity  1  or  2  has  Sch(U) 
=  2,  and  with  CNOT  complexity  0  has  Sch(U)  =  1 . 

Having  considered  CNOT  complexity  in  relation  to  Schmidt  number  for  two-qubit 
operators,  an  alternative  way  to  express  our  classification  is  in  terms  of  the  Hartley 
strength,  KHar(Q)  =  log2[Sch(Q)].  We  find  that  KHar(U)  =  2  for  operators  U  of  CNOT 
complexity  2  or  3,  KHar(U)  =  1  for  those  operators  with  CNOT  complexity  1  or  2,  and 
finally  KHar(U)  =  0  for  local  unitaries.  Among  examples,  the  particular  swap  operation  Us 
and  CNOT  itself  have  KHar=  1-  The  SWAP”  gate,  with  0  <  a  <  1,  as  well  as  F,  the 
quantum  Fourier  transform  on  two  qubits,  have  KHar=  2  and  CNOT  complexity  3.  The 
classification  that  we  have  given  is  hierarchical,  in  the  sense  that  operators  of  higher 
CNOT  complexity  or  Schmidt  number  may  simulate  those  of  lower  complexity,  but  not 
the  other  way  around. 

We  have  discussed  alternative  means  for  obtaining  both  the  CNOT  complexity  of  a  two- 
qubit  operator  and  its  canonical  decomposition  [coffeydeiotte,  cts].  Besides  the 
operational  decompositions  given  elsewhere,  we  have  illustrated  the  constructive 
procedure  of  Childs  et  al.  [childs]  for  finding  the  three  nonlocal  parameters  of  a  canonical 
decomposition. 
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The  operator  Schmidt  decomposition  has  with  it  a  normalization  that  is  convenient  for 
probability  and  entropy  considerations.  For  unitaries  acting  on  systems  A  and  B  of 
dimensions  dA  and  ds,  respectively,  the  relation  tr(U*U)=dAdB  gives  that  the  Schmidt 
coefficients  sj  satisfy  Sj  sj^  =  dAds.  Thus,  the  normalized  coefficients  Sj^/dAds  form  a 
probability  distribution.  For  the  two-qubit  operators  U  considered  here,  simply  the 
numbers  Sj  /4  give  a  probability  distribution.  The  Schmidt  strength  may  be  defined  as  the 
Shannon  entropy  of  this  distribution,  providing  a  measure  of  the  nonlocal  content  of  U. 

3,4  Scheduling  and  logistics  problems  via  adaptive  quantum  search 

The  recent  developments  in  the  field  of  quantum  computing  have  allowed  computer 
scientists  and  others  to  view  optimization  problems  from  a  new  perspective.  By  framing 
optimization  as  a  problem  of  unordered  database  search,  these  problems  can  be  solved 
using  algorithms  based  on  Grover's  quantum  search,  theoretically  providing  quadratic 
speedup  in  runtime.  In  1996  Grover  [grover,  groverl,  grover2,  groverS]  used  the 
property  of  quantum  parallelism  to  design  an  algorithm  to  search  an  unordered  database 
in  0(a/N)  time,  a  problem  that  classically  takes  0(N)  time  where  N  is  the  database  size. 
It  has  been  hypothesized  that  Grover's  results  can  be  generalized  to  implement  adaptive 
search,  an  algorithm  useful  to  optimization  which  cannot  be  implemented  efficiently  on 
classical  machines.  Combining  Grover  adaptive  search  with  quantum  encoding 
techniques,  it  may  be  possible  to  provide  better  than  quadratic  speedup  in  optimizing 
some  families  of  scheduling  problems. 

Scheduling  algorithms  attempt  to  solve  “scheduling  problems”,  which  consist  of  finding 
the  optimal  order  and  distribution  of  a  set  of  tasks  on  a  set  of  machines  or  other  resources. 
Classically,  this  class  of  problem  is  NP-hard,  meaning  that  it  cannot  be  solved  in 
polynomial  time,  and  practical  algorithms  to  the  problem  set  typically  provide  a  “good” 
solution  rather  than  the  optimal  solution  [lu].  Application  of  Grover's  quantum  search 
[grover]  to  these  problems  typically  still  cannot  produce  polynomial-time  algorithms; 
however,  it  can  provide  improved  exponential-time  algorithms  that  are  more  efficient 
than  classical  solutions.  For  example,  the  algorithm  proposed  in  [lu]  claims  to  return  a 
schedule  for  running  N  jobs  on  M  non-homogeneous  machines  as  well  as  the  longest  job 
runtime  in  0(a/M^)  time.  This  problem  is  known  as  the  R||Cmax  problem,  and  classically 
requires  0(M’^)  time  to  solve.  Lu  and  Marinescu  describe  in  [lu]  a  systematic  approach  to 
reformulating  other  scheduling  problems  so  that  they  can  take  advantage  of  Grover's 
search. 
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Scheduling  problems  emerge  in  a  wide  variety  of  scenarios.  One  example  is  in  operating 
system  design,  where  computer  resources  must  be  allocated  for  multiple  jobs  on  a  system. 
Another  example  is  dividing  tasks  among  workers  in  a  business.  Employees  have  varied 
abilities,  and  some  are  suited  better  to  a  certain  type  of  job  than  others.  Scheduling 
algorithms  ean  help  managers  decide  who  should  be  working  on  what  and  when.  The 
transportation  industry  faces  a  variety  of  seheduling  problems,  which  include  the 
additional  complexity  of  requiring  airplanes,  trains,  and  trucks  to  follow  connected  routes. 
As  quantum  computers  become  commercially  viable,  quantum  scheduling  algorithms  can 
be  used  to  solve  practical  scheduling  problems  which  are  very  difficult  and  time 
consuming  to  solve  classically. 

Before  discussing  adaptive  search,  we  describe  some  necessary  background  on  Grover's 
search  algorithm  [grover,  groverl,  grover2,  kaye,  nielsen].  Grover's  search  is  an 
algorithm  to  search  an  unordered  database.  This  problem  is  like  trying  to  search  a  phone 
book  to  find  which  name  corresponds  to  a  phone  number.  Since  there  is  no  structural 
information  to  narrow  down  the  search,  classically  each  entry  in  the  phone  book  must  be 
ehecked  sequentially  to  see  if  the  assoeiated  entry  matches  the  queried  number.  In  the 
average  case,  the  algorithm  will  search  half  of  the  entries  before  finding  the  target.  In  the 
worst  case,  there  is  no  target  (the  number  belongs  to  someone  not  in  the  phone  book),  and 
the  algorithm  searches  the  whole  database  before  returning.  Grover's  algorithm  takes 
advantage  of  a  quantum  computer's  ability  to  manipulate  a  superposition  of  input  states 
simultaneously  and  can  find  the  target  entry  in  0(a/N)  “lookups”,  where  N  is  the  number 
of  elements  in  the  database.  The  notion  of  a  lookup  (oracle  funetion  query)  is  slightly 
different  between  the  quantum  mechanieal  and  the  classical  techniques  since  the  quantum 
algorithm  can  “check”  all  of  the  input  states  at  once.  In  some  cases,  there  may  be 
multiple  target  states  in  the  database,  and  a  suecessful  search  would  return  any  of  these 
states. 

The  first  step  of  the  algorithm  is  to  initialize  the  system  to  an  equal  superposition  of  all 
input  states.  Typically,  N  is  assumed  to  be  a  power  of  2,  so  N  =  2”.  Since  memory  is 
usually  available  in  powers  of  2,  this  assumption  simplifies  the  algorithm  because  it  does 
not  need  to  account  for  unused  states. 
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Next,  the  “Grover  gate”  is  applied  [7t:a/N}/4]  times.  The  Grover  gate  is  made  of  two  main 
eomponents,  a  seleetive  phase  shift  Uf  and  an  “inversion  about  average”  funetion  Uv,,. 
The  seleetive  phase  shift  rotates  the  target  states  by  n,  essentially  inverting  their 
probability  amplitudes.  This  funetion  is  known  as  the  “oraele”  funetion,  whieh  is  the 
quantum  equivalent  of  the  “lookup”  funetion.  Unlike  the  elassieal  version,  it  does  not 
direetly  indieate  what  a  target  value  is,  but  merely  affeets  the  probability  amplitude  of  the 
target  states.  The  inversion  about  average  operation  refleets  the  probability  amplitude  of 
eaeh  state  about  the  average  amplitude.  These  two  operations  inerease  the  amplitude  of 
the  target  state  and  reduee  the  probability  of  all  other  states.  Sinee  the  probability  of 
measuring  a  partieular  value  is  the  eorresponding  state's  amplitude  squared,  when  a 
measurement  is  performed,  a  target  value  will  be  observed  with  probability  greater  than 
1/2.  Grover's  seareh  does  not  return  a  target  state  with  probability  1,  but  if  it  is  performed 
more  than  onee,  the  probability  that  the  result  returned  is  a  target  state  is  inereased.  The 
returned  result  eould  also  be  eheeked  elassieally  to  verify  eorreetness,  although  this 
eonfirmation  is  not  applieable  in  all  oases. 

While  the  number  of  Grover  iterations  required  to  maximize  the  amplitude  of  the  target 
state(s)  is  less  than  0(a/N),  the  amplitude  fluotuates  signifioantly  based  on  the  number  of 
iterations  and  the  proportion  of  the  states  that  are  marked.  The  probability  of  finding  a 
target  state  is  given  by: 

gr(p)  =  sin^[(2r  +  l)arosin(A/p)]  (5) 

where  p  is  the  fraotion  of  states  that  are  marked  and  r  is  the  number  of  applioations  of  the 
Grover  gate  [baritompa].  A  plot  of  g6(p)  is  shown  in  Figure  3.  Depending  on  the  value 
of  p,  6  applioations  of  the  Grover  gate  eould  be  suooessful  with  either  very  high,  or  very 
low  probability.  If  the  ratio  p  was  known,  one  eould  easily  determine  the  best  number  of 
iterations  to  apply,  but  p  is  typioally  unknown. 
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Figure  3.  The  probability  distribution  of  observing  a  target  state  after  6  applications  of  the  Grover 
gate  with  the  fraction  p  of  target  states.  Depending  on  the  value  of  p,  the  search  success  probability 

varies  greatly. 


3.4.1  Adaptive  Search 

Adaptive  search  is  an  algorithm  for  finding  the  minimum  value  of  f(x)  for  a  finite  set  of 
values  X  e  X  by  selecting  decreasing  values  f(x)  until  the  minimum  is  reached.  The 
algorithm  is  as  follows  [bulger]; 

X  =  random  element  of  X 

y  =  f(x) 

FOR  {i  =  1,2...  while  f(xi)  is  not  the  minimum  of  f(x)  for  x  e  X} 

Xi  =  {x  e  X  I  f(x)  <  Yi} 
x'  =  random  element  of  Xi 
y’  =  f(x’) 

Xi+i  =  x'  and  yi+i  =  y' 

ENDFOR 
return  x' 
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There  are  two  major  problems  with  implementing  this  algorithm  classieally.  First,  the 
terminating  condition  is  that  f(xi)  is  not  the  minimum,  meaning  that  in  order  to  find  value, 
X,  that  minimizes  f(x),  the  minimum  of  f(x)  must  already  be  known.  Second,  identifying 
the  set  Xi,  known  as  the  improving  region,  classically  requires  checking  each  element  of 
the  previous  improving  region,  negating  any  potential  speedup  of  adaptive  search. 
Grover's  search,  however,  provides  remedies  for  both  of  these  flaws. 

A  quantum  algorithm  to  find  the  minimum  of  a  finite  set  was  first  proposed  by  Durr  and 
Hoyer  in  [durr].  They  proposed  replacing  the  stopping  condition  by  a  time  limit  restraint. 
According  to  their  analysis,  after  running  for  22.5a/N  +  1.41og^N  units  of  time,  their 
algorithm  would  return  the  minimum  value  with  probability  of  greater  than  1/2. 
Baritompa  et  al.  determined  that  this  result  had  incorrect  constants,  but  was  correct  in 
order  of  magnitude,  providing  more  rigorous  analysis  considering  other  parameters  of  the 
search  [baritompa]. 

The  second  problem  is  solved  naturally  using  the  Grover  operator.  Using  the  oracle 
function,  h(w)  =  (f(w)  <  Y),  the  target  states  are  exactly  the  states  in  the  improving  region. 
A  successful  application  of  Grover  search  with  this  oracle  function  will  return  a  random 
element  from  the  improving  region  with  uniform  distribution.  As  Xi  becomes 
increasingly  close  to  the  optimal  value,  however,  the  improving  region  shrinks  and  more 
Grover  iterations  are  potentially  required  for  the  search  to  be  successful. 

Grover's  search  introduces  a  new  complication  to  adaptive  search  due  to  the  probability 
distribution  defined  by  gr(p).  Each  successful  iteration  of  adaptive  search  will  decrease 
the  fraction  p,  having  a  potentially  dramatic  effect  on  the  probability  of  search  success.  If 
the  number  of  applications  of  the  Grover  gate  does  not  change  between  iterations,  it  is 
possible  to  stall  at  a  point  of  low  success  probability,  the  minima  of  Figure  3.  Varying 
the  value  of  r  can  potentially  avoid  this  problem  and  several  methods  of  “rotation 
scheduling”  have  been  proposed.  The  first  is  to  select  a  random  value  between  0  and  m 
where  m  is  initially  1,  and  is  multiplied  by  a  factor  A,  after  each  failed  search  attempt 
[baritompa].  This  approach  has  the  advantage  of  being  simple  to  compute,  but  m  grows 
exponentially  with  the  number  of  failures.  Baritompa  et  al.  propose  a  predefined 
sequence  of  rotations  based  on  an  estimate  of  p  after  each  iteration.  This  sequence  of  r 
values  performs  slightly  better  in  some  cases,  but  computing  the  sequence  requires 
exponential  time  and  this  sequence  performs  equally  well  with  the  randomized  sequence 
when  there  are  repeated  values  in  the  range  in  the  search  domain.  Neither  approach 
guarantees  search  success  as  p  approaches  1/N. 

Adaptive  search  has  several  important  properties  for  optimization.  First,  unlike  many 
classical  minimization  algorithms,  adaptive  search  is  not  affected  by  local  minima.  Since 
the  improving  region  is  defined  by  the  range  of  the  function,  not  the  domain,  the  global 
minimum  will  never  be  excluded  from  the  improving  region.  Second,  while  most  of  the 
work  done  by  the  Grover  implementation  of  adaptive  search  is  done  during  the  final 
iteration  when  the  improving  region  is  as  small  as  possible,  the  intermediate  results  are 
generally  good  approximations  of  the  minimum,  especially  if  there  are  multiple  values 
only  slightly  greater  than  the  global  minimum. 
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3.4.2  Framing  Optimization  Problems  as  Search  Problems 
The  concept  of  quantum  parallelism  makes  it  possible  to  treat  an  optimization  problem  as 
a  search  problem  with  far  fewer  resources  than  are  classically  required.  Classically,  a 
brute  force  approach  to  optimization  would  be  to  test  the  set  of  all  possible  outcomes  of 
the  system  being  optimized  (the  state  space  of  the  system).  Since  the  state  space  typically 
grows  exponentially  in  size  as  systems  become  more  complex,  searching  the  entirety  of 
the  state  space  is  usually  computationally  intractable.  To  overcome  this,  pruning  is 
typically  applied  to  reduce  the  size  of  the  state  space  by  excluding  clearly  non  optimal  or 
duplicate  results.  Correct  and  effective  pruning  requires  knowledge  of  attributes  and 
symmetries  of  the  system  that  are  often  either  unknown,  or  may  not  exist. 

Classically,  each  state  in  the  state  space  must  be  computed  separately,  either  sequentially, 
or  with  a  low  level  of  parallellization.  In  a  quantum  computer,  however,  the  entire 
statespace  can  be  computed  in  a  superposition  with  one  operation.  Once  the  state  space 
has  been  prepared,  algorithms  like  Grover  adaptive  search  can  be  applied  to  determine 
the  optimum  solution. 


Lu  and  Marinescu  propose  in  [lu]  a  procedure  for  generating  the  state  space  for  the  R  || 
Cmax  problem,  which  appears  to  be  easily  adaptable  to  other  optimization  situations.  The 
quantum  circuit  described  makes  use  of  the  entanglement  between  qubits  by  using 
controlled  gates  triggered  by  a  set  of  index  qubits.  A  circuit  computing  the  state  |i)|Gi)  is 
shown  in  Figure  4.  If  the  input  to  the  index  qubits  is  an  equal  superposition  of  states, 
which  can  easily  be  prepared  with  Hadamard  gates,  and  if  all  the  necessary  gates  can  be 
implemented  efficiently,  circuits  like  this  can  be  chained  together  to  create  an  equal 
superposition  of  all  states  in  the  state  space. 


|i>  (g)  |G|> 


Figure  4.  A  quantum  circuit  computing  the  state  |i)|Gi).  The  index  lines  ensure  that  the  gates  Gi  are 

each  applied  only  once. 
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3.4.3  Using  Grover  Search  to  Find  the  Minimum 

Once  the  state  space  has  been  computed,  Grover  adaptive  search  can  be  used  to  find  the 
minimum  or  maximum  of  a  particular  parameter  encoded  in  the  superposition.  Using  an 
oracle  function  h(w)  defined  as  h(w)  =  (f(w)  <  Y)  where  Y  is  an  adjustable  threshold 
value  and  w  is  the  set  of  qubits  of  interest,  an  adaptive  search  like  the  one  described  in 
[baritompa]  can  be  used  to  find  the  optimal  value.  The  encoding  described  by  Lu  and 
Marinescu  entangles  the  input  states  with  the  output  states,  so  when  a  measurement  is 
finally  performed,  the  input  states  collapse  into  the  set  of  inputs  producing  the  measured 
output.  This  algorithm  returns  not  only  the  optimal  output  value,  but  how  to  achieve  it. 
Using  a  different  encoding  circuit  or  oracle  function  would  produce  a  circuit  that  could 
minimize  or  maximize  a  variety  of  parameters  of  a  system. 

3.4.4  Results:  Adaptive  search  simulation 

In  order  to  analyze  the  properties  of  adaptive  search,  several  simulations  were  developed 
using  classical  techniques.  These  simulations  verified  the  desirable  properties  of  adaptive 
search.  For  example,  the  ID  function  shown  in  Figure  5  demonstrates  that  the  algorithm 
is  not  affected  by  local  minima.  Also,  after  a  small  number  of  iterations,  all  subsequent 
values  are  within  a  small  range  of  the  optimal  value,  even  though  they  are  in 
disconnected  parts  of  the  domain. 


Figure  5.  Successive  iterations  of  adaptive  search  (circles)  on  a  ID  fnnction  (solid  line)  with  mnltiple 
local  minima.  The  search  converges  qnickly  to  the  optimnm  valne  despite  local  extrema. 
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Extending  this  simulation  to  a  2D  function  already  begins  to  demonstrate  the  limits  of 
classical  implementations  of  adaptive  search.  Since  the  domain  grows  exponentially  in 
the  number  of  dimensions,  each  additional  dimension  requires  exponentially  more 
resources  to  find  the  optimal  value.  However,  in  2  dimensions,  the  same  properties  hold, 
as  shown  in  Figure  6  and  Figure  7. 


Figure  6.  A  3D  contour  plot  of  a  2D  function  with  successive  adaptive  search  iterates.  2D  adaptive 
search  exhibits  the  same  desirable  properties  as  the  ID  variant. 

4 

xIO 


7- 


5- 
4- 
3- 
Z- 
1  ^ 


Figure  7.  Adaptive  search  applied  to  the  Rosenbrock  banana  function,  a  standard  optimization 

objective  function  [bananaFunc]. 
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3.4.5  Rotation  Schedule  Implementation 


Baritompa  et  al.  provide  in  [baritompa]  an  algorithm  for  eomputing  a  sequenee  of 
rotation  counts  to  improve  the  probability  of  suecess  of  sequential  adaptive  search 
iterations.  The  algorithm  requires  indefinite  integration  of  successive,  increasingly 
eomplieated  funetions,  whieh  are  not  easily  integrated.  One  way  to  aetually  eompute 
these  integrals  is  to  represent  the  funetion  gr(p)  as  a  polynomial.  The  algorithm  to 
generate  the  rotation  sehedule  was  implemented  using  a  Taylor  series  approximation  to 
gr(p)  as  well  as  an  exaet  polynomial  representation  based  on  the  Gauss  hypergeometrie 
funetion  2F1, 


/ 


gXp)  =  {'2-r  +  lf 


r  +  \-r\-\p 
V  2  y 


(6) 


Examples  are  gi(p)  =  (3  -  4p)^p  =  9p  -  24p^  +  16p^  and  g2(p)  =  p[5  +  4p(4p  -  5)]^.  The 
degree  of  the  terminating  2F 1  function  in  Equation  6  is  r,  so  that  the  overall  degree  of  gr(p) 
in  p  is  2r  +  1.  Due  to  the  oseillatory  nature  of  gr(p),  shown  in  Figure  3,  a  very  high  order 
Taylor  polynomial  is  required  as  r  grows,  rendering  this  approaeh  infeasible  at 
suffieiently  large  r.  The  hypergeometrie  funetion-based  approaeh  does  not  suffer  from 
this  problem,  but  as  the  integrals  become  inereasingly  eomplex,  the  computation  time 
inereases  dramatieally.  With  computational  resources  readily  available  (desktop  PC),  the 
hypergeometrie  funetion  implementation  was  able  to  eompute  the  first  30  terms  of  the 
rotation  sehedule  in  24  hours  in  Mathematiea.  A  plot  of  the  first  30  terms  of  the  sequenee 
is  shown  in  Figure  8. 


r  n 


Figure  8.  The  first  30  terms  of  the  rotation  schedule  descrihed  hy  [haritompa]. 
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3.4.6 


Future  Work 


3.4.6.1  Adaptive  Search 

The  usefulness  of  Grover  adaptive  seareh  hinges  on  the  hypothesis  that  an  oraele  funetion 
like  h(w)  =  (f(w)  <  Y)  ean  be  implemented  effieiently.  More  researeh  is  required  to 
design  and  effieiently  implement  a  quantum  eireuit  that  ean  implement  eomparison 
operators  such  as  “greater  than”  or  “less  than”.'  This  gate  must  be  able  to  take  two  inputs 
because  the  threshold  value  Y  will  change  with  each  iteration,  so  a  hard  coded  value  will 
not  suffice.  It  remains  an  open  question  just  how  efficiently  the  needed  oracles  for 
generic  problems  may  be  implemented. 

Baritompa  et  al.  analyzed  the  Grover  adaptive  search  algorithm  to  determine  the  expected 
number  of  iterations  before  finding  the  optimal  value  [baritompa].  However,  a  bound  on 
how  close  to  the  optimal  value  can  be  expected  after  a  number  of  iterations  would  be 
valuable  as  well.  For  example,  if  the  second  most  optimal  value  can  be  found  in  half  the 
time  required  to  find  the  optimal  value,  it  may  be  unnecessary  to  perform  as  many 
iterations. 

3.4.6.2  Scheduling  Problem  Encoding 

The  encoding  system  developed  in  [lu]  appears  to  be  valid,  however,  there  are  a  number 
of  extensions  that  would  likely  be  required  in  practice.  The  circuit  as  is  requires  that  all 
inputs  be  integers.  This  is  not  a  limitation  in  practice,  since  the  least  common  multiple  of 
all  necessarily  rational  input  parameters  may  be  used. 

3.4.6.3  Termination  of  the  algorithm 

Setting  a  run  time  bound  on  the  size  of  the  improving  region  to  the  total  domain  is  a 
working  possibility  for  a  criterion  for  terminating  the  algorithm.  However,  many 
common  global  minima  could  be  problematic.  We  may  note  that  even  in  the  easier  “find 
Min”  situation  of  Durr  and  Hoyer  there  is  an  interrupt  on  run  time.  The  main  question  is 
again  how  well  this  can  be  done  for  generic  problems. 

3.4.6.4  Combining  encoding  and  search 

Once  the  quantum  encoding  of  a  scheduling  problem  is  shown  to  be  correct,  and  adaptive 
search  has  been  shown  to  be  as  efficient  as  predicted,  the  two  must  be  combined.  In 
order  to  show  that  scheduling  problems  can  actually  be  solved  using  this  technique,  the 
whole  system  must  be  simulated  and  benchmarked.  The  number  of  quantum  gates 
required  to  construct  these  circuits  could  also  be  intractably  large,  so  optimization  of  the 
circuit  implementation  may  also  be  necessary  to  realize  this  algorithm. 
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3.4.6.5  Concluding  thoughts  on  adaptive  search 

As  quantum  computers  become  closer  to  being  viable,  the  development  of  quantum 
algorithms  is  becoming  increasingly  important.  A  quantum  computer  has  the  potential  to 
solve  problems  such  as  scheduling  problems  that  elassically  require  prohibitive  resources. 
Solutions  to  scheduling  problems  would  have  a  broad  reaching,  positive  impact  on  a 
variety  of  industries.  Techniques  exist  to  solve  pieces  of  scheduling  problems  on  a 
quantum  computer:  a  quantum  encoding  of  the  problem,  and  an  algorithm  for 
optimization.  By  combining  these  teehniques,  it  may  be  possible  to  solve  scheduling 
problems  that  are  intraetable  using  elassieal  techniques. 


3.5  Feynman  diagram  integrals  and  special  functions 

Perturbative  quantum  field  theory  visually  and  systematieally  describes  higher  order 
corrections  by  means  of  Feynman  diagrams.  In  many  cases  the  corresponding 
multidimensional  Feynman  integrals  may  be  reduced  all  the  way  to  one  dimension.  In 
recent  years  both  rigorously  and  empirically  based  expressions  for  these  integrals  in 
terms  of  values  of  special  functions  have  appeared.  In  a  series  of  articles,  we  have  made 
a  number  of  contributions  for  closed  form  results  for  redueed  Feynman  integrals 
[coljmpl,  coljmp2,  cj,  cl,  coffeyOS].  We  were  also  able  to  advance  the  state  of  the  art  of 
the  theory  for  particular  special  functions  of  mathematical  physics,  including  the  Clausen 
function  CI2  and  the  dilogarithm  function  Li2.  Other  speeific  speeial  functions  that  were 
very  useful  included  the  generalized  hypergeometric  functions  and  the  polylogarithms. 
Results  to  date  indicate  connections  between  multiple  subjects  including  hyperbolic 
geometry,  analytic  number  theory,  special  functions,  and  Feynman  diagrams.  A  very 
recent  example  [coljmpl]  evaluates  a  highly  symmetrie  Feynman  integral  in  terms  of  a 
remarkably  compact  difference  of  Clausen  function  CI2  values,  proving  a  numerically 
based  conjecture  that  was  open  for  approximately  ten  years. 

A  certain  dilogarithmic  integral  I7  turns  up  in  a  number  of  eontexts  ineluding  Feynman 
diagram  calculations,  volumes  of  tetrahedra  in  hyperbolic  geometry,  knot  theory,  and 
conjectured  relations  in  analytic  number  theory.  We  provided  an  alternative  explicit 
evaluation  of  a  parameterized  family  of  integrals  eontaining  this  partieular  case.  By 
invoking  the  Bloch-Wigner  form  of  the  dilogarithm  function,  we  produced  an  equivalent 
result,  giving  a  third  evaluation  of  I7.  We  also  alternatively  formulated  some  conjectures 
which  we  posed  in  terms  of  values  of  the  specific  Clausen  function  CI2  [coljmp2, 
coffeyOS]. 
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The  evaluation  of  harmonic  number  sums  has  been  useful  in  several  areas  of  mathematics, 
including  analytic  number  theory,  for  some  time.  More  recently,  the  evaluation  of  Euler 
sums  has  been  shown  important  in  various  areas  of  theoretical  physics,  including  in 
support  of  Feynman  diagram  calculations.  Even  more  recently,  it  has  been  shown  that 
the  evaluation  of  generalized  harmonic  number  sums  is  very  useful  in  resolving  open 
questions  on  Feynman  diagram  contributions  and  relations  among  special  functions 
[coljmpl],  including  the  dilogarithm,  Clausen  function,  and  generalized  hypergeometric 
function.  Harmonic  number  sums  also  occur  in  computer  science  in  the  efficiency 
analysis  of  algorithms.  Especially  in  the  analysis  of  sorting  and  searching  algorithms, 
harmonic  number  sum  evaluation  and  asymptotic  analyses  are  useful. 

Motivated  by  such  applications,  we  have  developed  results  on  generalized  harmonic 
number  sums  [cl].  We  give  results  where  evaluations  are  possible  in  terms  of  generic 
polylogarithm  functions  Elk,  although  we  give  specializations  to  the  dilogarithmic  and 
trilogarithmic  cases,  where  a  fuller  body  of  theory  exists,  and  relatively  more  instances  of 
explicit  expression  of  specific  values  in  terms  of  elementary  functions  is  possible.  We 
demonstrated  the  evaluation  of  a  class  of  generalized  harmonic  number  sums,  presented 
examples,  and  then  gave  representative  results  for  obtaining  identities  among  harmonic 
number  sums  in  terms  of  polylogarithmic  functions. 

As  a  result  of  these  research  activities,  the  PI  was  invited  as  a  speaker  at  a  very  recent 
workshop  dealing  with  Feynman  diagrams,  multiple  zeta  values,  and  integrals 
[workshop].  Other  areas  of  mathematics  covered  in  this  workshop  were  abstract  algebra, 
combinatorics,  and  graph  theory.  The  main  areas  of  theoretical  physics  concerned  were 
quantum  field  theory  and  quantum  statistical  mechanics.  The  PI  gave  two  lectures  at  the 
graduate  level  on  Feynman  diagrams,  integration,  and  special  functions. 

3,6  QLGAs  for  the  Maxwell  equations 

In  our  research  we  showed  that  a  quantum  lattice  gas  approach  can  provide  a  viable 
means  for  numerically  solving  the  classical  Maxwell  equations.  By  casting  the  Maxwell 
equations  in  Dirac  form,  the  propagator  may  be  discretized,  and  we  described  [coffey] 
how  the  accuracy  relative  to  the  time  step  may  be  systematically  increased.  The  quantum 
lattice  gas  form  of  the  discretization  is  suitable  for  implementation  on  hybrid  classical- 
quantum  computers.  In  the  published  paper  [coffey]  we  also  discuss  a  number  of 
extensions,  including  application  to  inhomogeneous  media. 

The  analytic  and  numerical  modeling  of  electromagnetic  phenomena  is  a  subject  of  much 
practical  importance.  Due  to  a  panoply  of  ingredients  including  geometry, 
dimensionality,  boundary,  initial,  and  radiation  conditions,  and  constitutive  relations, 
nature  presents  a  rich  variety  of  these  phenomena. 
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In  this  section  we  introduce  the  formulation  of  the  Maxwell  equations  as  quantum  lattice 
gas  algorithms.  In  this  way,  the  partial  differential  equations  of  electromagnetism  are 
discretized  in  both  space  and  time,  permitting  numerical  simulation.  Here  we  concentrate 
on  the  theoretical  framework,  but  we  provide  sufficient  detail  that  software 
implementations  could  be  developed.  Many  extensions  of  our  approach  are  possible,  and 
we  described  several  of  them  elsewhere  [coffey]. 

We  exploit  the  formulation  of  the  Maxwell  equations  for  the  electromagnetic  field  in 
terms  of  a  multicomponent  wave  function,  and  this  wave  function  need  not  correspond  to 
a  spin  1  particle.  Auxiliary  conditions  are  used  as  necessary  to  ensure  that  the  full  set  of 
Maxwell  equations  is  satisfied. 

It  has  been  known  for  some  time  that  the  Maxwell  equations  may  be  written  in  Dirac 
form  [iwoprd,  moses].  We  complemented  this  framework,  and  provided  supporting 
techniques  so  that  competitive  algorithms  could  be  developed.  Generally  the  Maxwell 
equations  in  the  continuous  setting  may  be  brought  into  the  form: 

5t\|/  =  K  S  •V\|/,  (7) 

where  k  is  a  constant,  Sj  are  three  spin  matrices  of  suitable  dimension,  and  \\f  the 
corresponding  wave  function.  We  focused  on  V  as  the  three  dimensional  (3D)  gradient 
operator,  but  this  is  by  no  means  necessary.  The  result  of  discretizing  Equation  7  is  to 
present  the  time  evolution  of  \)/  in  terms  of  spatial  shift  and  spin-component  mixing 
operations.  In  the  language  used  in  a  very  related  setting,  these  operations  may  be 
rewritten  as  streaming  and  collision  operations,  respectively. 

Having  written  a  discretization  of  Equation  7,  we  take  the  further  step  of  writing  the 
discrete  evolution  operator  in  terms  of  a  product  of  collision  and  streaming  operators.  It 
is  this  form  of  the  discrete  evolution  that  could  be  implemented  on  a  hybrid  processor, 
wherein  local  nodes  of  a  few  qubits  are  connected  via  classical  communication  channels 
to  neighboring  nodes.  The  on-site  collision  operator  generates  entanglement  within  the 
individual  nodes.  As  it  is  unitary,  it  can  be  realized  in  a  quantum  system,  and  further 
provides  a  stable  algorithm.  In  sum,  our  procedure  shows  how  to  discretize  systems  such 
as  Equation  7  in  a  number  of  settings,  and  importantly  in  a  factorized  form  suitable  to 
combined  classical-quantum  computing. 
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Lattice  gas  algorithms  have  proved  very  effieient  and  versatile  in  a  variety  of  eontexts, 
ineluding  fluid  dynamies  and  plasma  physios  and  other  multi-partiole  simulations 
[rothman,  yepez,  vahala].  They  have  been  used  to  numerioally  solve  the  Navier  Stokes 
and  Boltzmann  equations  in  diverse  geometries,  upon  a  range  of  oomputing  platforms.  A 
key  feature  of  these  algorithms  is  their  suitability  for  parallel  and  distributed  oomputing 
(e.g.,  [harting]).  We  point  out  that  it  should  be  possible  to  develop  these  same 
advantages  for  solving  the  Maxwell  equations.  Moreover,  there  are  now  prospeots  for 
nearer  term  hybrid  olassioal-quantum  oomputing  (e.g.,  [pravia,  hems,  ohen06])  as  well  as 
more  distant  purely  quantum  prooessors  [nielsen].  Various  quantum  lattioe  gas 
algorithms  have  been  proposed  for  the  solution  of  the  Dirao  and  Sohroedinger  equations 
[boghosian98,  yepez02].  Further  development  of  the  types  of  algorithms  that  we  suggest 
would  allow  the  exploitation  of  entanglement  in  quantum  hardware.  This  would  be  an 
additional  resouroe,  on  top  of  the  olassioal  parallelization  aspeets  of  these  algorithms. 

Signifieant  portions  of  our  work  made  use  of  results  of  I.  Bialynieki-Birula  [iwoprd, 
iwoaeta,  iwohep],  who  was  mainly  interested  in  theoretieal  eonstruetions.  In  [eoffey]  we 
reealled  the  formulation  of  the  Maxwell  equations  in  Dirao  (spin  1/2)  form,  providing  a 
self-0 ontained  development,  whioh  may  make  this  subjeot  aooessible  to  a  wider  audienoe. 
We  then  oonsidered  a  spin  1  formulation  of  the  Maxwell  equations  and  again  took  it  as 
the  basis  of  a  oomputational  method.  For  both  formulations,  we  provided  supporting 
details.  We  then  desoribed  a  number  of  extensions,  ineluding  how  to  improve  the 
oonvergenoe  rate  of  the  disoretizations,  this  latter  important  topio  applies  to  either  the 
spin  1/2  or  spin  1  formulations. 

We  have  shown  how  the  multieomponent  wave  equations  of  relativistie  quantum 
meehanies  may  be  adapted  for  a  deseription  and  numerieal  solution  of  the  elassieal 
Maxwell  equations.  Written  in  a  spin  1  form,  the  wave  funetion  has  three  independent 
eomponents  that  are  direetly  applieable  to  deseribing  the  massless  photon,  whereas 
written  in  a  spin  1/2  Dirae  form,  the  four-eomponent  wave  funetion  requires  a  eonstraint 
eondition.  By  diseretizing  the  quantum  evolution  operator  (propagator),  we  have  derived 
algorithms  for  integrating  the  Maxwell  equations.  These  algorithms  provide  an  explieit 
numerieal  integrator  and  are  suitable  to  any  number  of  spatial  dimensions.  This  approaeh 
eontrasts  with  that  of  first  diseretizing  the  system  Hamiltonian,  whioh  may  often  lead  to 
physioal  and  mathematioal  diffioulties. 

As  seen  in  detail  in  [eoffey],  there  are  equivalenoes  with  effeotive  finite  differenoe 
approximations  of  the  partial  differential  equations,  and  this  is  not  surprising 
mathematioally.  However,  physioally  our  approaeh  provides  oonneotions  with  the  path 
integral  formulation  of  quantum  meehanies,  and  this  was  a  motivating  reason  behind  the 
work  of  [iwoprd].  Indeed,  an  iteration  of  a  disoretized  time  evolution  gives  a  sum  over 
Dirao  partiole  histories. 
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We  emphasize  that  mathematieal  equivalenees  are  quite  distinet  from  software  and 
hardware  implementations.  Whereas,  to  be  very  speeifie,  a  finite  differenee 
representation  is  suitable  for  either  serial  or  parallel  proeessing,  the  lattiee  gas  forms  are 
also  very  suited  to  a  parallel  arehitecture.  The  lattice  gas  representation  comprises  two 
main  stages.  One  is  the  local  aspect,  embodied  in  W  or  collision  matrices  C  acting  on 
local  field  or  wave  function  components.  The  other  is  the  advective  (or  streaming)  stage 
shown  by  the  shift  operators  moving  field  information  to  neighboring  sites.  Collections 
of  lattice  sites  may  be  mapped  to  local  computing  nodes  of  minimal  resources  within  a 
parallel  machine.  For  instance,  for  a  square  lattice,  blocks  of  rows  or  columns  of  the 
lattice  may  be  mapped  to  individual  processors,  or  instead  checkerboard  subblocks  may 
be  mapped  to  the  processors.  In  this  way,  parallel  computers  with  a  large  number  of 
nodes  but  each  with  limited  computing  capability  and  memory  may  be  taken  advantage 
of,  with  many  variations  on  these  ideas  possible. 

There  is  no  difficulty  in  increasing  the  rate  of  convergence  of  our  family  of  algorithms. 
One  has  only  to  apply  systematic  Trotter- Suzuki  formulas  [suzuki]  for  improved 
discretization  of  the  propagator. 

We  have  briefly  given  an  extension  of  the  lattice  gas  approach  to  inhomogeneous  media, 
and  expect  that  other  extensions  are  possible  [coffey].  Already  in  one  or  two  dimensions, 
this  opens  the  way  for  solving  a  variety  of  interesting  and  practical  problems. 

While  we  have  focused  on  the  Maxwell  equations,  many  other  applications  could  be 
developed  when  the  systems  of  governing  partial  differential  equations  can  be  cast  in  the 
form  of  Equation  7,  and  this  point  furthers  the  argument  for  the  very  wide  applicability  of 
lattice  gas  algorithms.  One  example  extended  coupled  system  is  the  Maxwell- 
Schroedinger  equations.  Here  the  absolute  square  of  the  wave  function  \|/  and  \|/*V\|/- 
v)/Vv)/*  feed  into  the  charge  and  current  densities,  respectively.  Another  example  coupled 
system  is  the  Maxwell-London  equations  for  a  phenomenological  description  of 
superconductivity,  wherein  the  London  relation  links  the  vector  potential  to  the 
supercurrent  density. 


3.7  Telegraph  Equation  QLGA 

The  telegraph  equation  combines  features  of  both  the  diffusion  and  wave  equations  and 
has  many  applications  to  heat  propagation,  transport  in  disordered  media,  image 
enhancement,  and  elsewhere.  In  this  section  we  give  an  overview  of  a  new  quantum 
lattice  gas  algorithm  for  this  partial  differential  equation  with  one  spatial  dimension.  This 
algorithm  generalizes  one  previously  known  for  the  diffusion  equation.  In  [cc2009]  we 
present  many  further  details,  including  an  analysis  of  the  algorithm  and  accompanying 
simulation  results.  The  QLGA  is  suitable  for  simulation  on  combined  classical-quantum 
computers. 
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Quantum  lattice  gas  algorithms  are  well  known  to  be  versatile  in  simulating  a  wide  range 
of  physieal  phenenomena.  Like  their  relatives  eellular  automata,  from  simple  loeal  rules, 
eomplex  dynamics  may  emerge  [brennen,  romanelli]. 

Lattiee  gas  algorithms  are  attraetive  due  to  their  relative  simplicity,  physical  foundations, 
and  suitability  for  implementation  on  parallel  eomputing  arehiteetures.  Lattiee  gas 
algorithms  may  ineorporate  eonservation  of  mass,  momentum,  and  energy,  and  in  the 
quantum  eontext,  probability.  Lattiee  gas  algorithms  have  proven  suceessful  in  a  range 
of  applieations  ineluding  fluid  dynamies  and  plasma  physics  and  other  multi-partiele 
simulations  [doolen,  rothman,  yepez], 

Reeent  experiments  and  proposals  for  eombined  elassical-quantum  eomputing  (e.g., 
[herns,  ehen06,  pravia])  further  motivate  the  development  of  quantum  lattice  gas  methods. 
For  instanee,  a  QLGA  for  the  linear  diffusion  equation  has  been  demonstrated  in  a  liquid- 
state  nuelear  magnetie  resonanee  system  [pravia].  In  addition,  a  detailed  design  for 
executing  a  QLGA  for  the  linear  diffusion  equation  with  superconducting  qubits  has  been 
given  [herns].  Sueh  implementations  eould  allow  the  exploitation  of  quantum 
entanglement  well  before  large-seale  purely  quantum  eomputers  are  eonstrueted. 

In  this  section  we  discuss  a  new  QLGA  for  simulation  of  the  telegraph  equation.  This 
hyperbolie  partial  differential  equation  combines  aspects  of  both  the  wave  and  diffusion 
equations.  As  sueh,  our  algorithm  subsumes  some  earlier  work  restrieted  to  the  diffusion 
equation.  Representative  numerieal  results  verifying  the  algorithm  and  its  analysis  are 
published  elsewhere  [ee2009]. 

Classieal  eonneetions  between  random  walk  and  the  telegraph  equation  have  been  known 
for  quite  some  time  [goldstein,  kac,  gaveau,  morette,  eodling].  In  sueh  a  model  on  a  one¬ 
dimensional  (ID)  lattiee,  a  walker  steps  a  distanee  Ax  in  time  inerement  t  randomly  to 
the  left  or  right,  with  additionally  a  probability  ar  to  reverse  direetion.  In  the 
simultaneous  limit  that  a  ^  oo  as  well  as  the  speed  v,  sueh  that  the  ratio  2a/v  =  1/D 
remains  eonstant,  the  diffusion  equation  results. 


du  _  ^  d^u 
dt  dx^ 


(8) 


Another  situation  for  whieh  this  equation  results  is  when  ar  =1/2.  Then  there  is  equal 
probability  for  a  move  to  the  left  or  right.  In  a  sense,  we  seek  a  quantum  version  of  this 
stoehastic  model. 
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We  also  mention  a  eonneetion  of  the  telegraph  equation  with  relativistic  quantum 
mechanics  and  the  point  of  view  of  a  Dirac  particle  as  moving  at  the  speed  of  light  c  with 
random  reversals  of  direction.  If  we  write  the  telegraph  equation  in  the  form  (5t  -2adrC 
dx^)P  =  0,  then  the  change  of  dependent  variable  P(x,t)  =  e'^  V(xd)  shows  that  \|/  satisfies 
the  Klein-Gordon  equation  (dt  -c  5x  -a  )\i/  =  0.  Both  the  telegraph  and  Klein-Gordon 
equations  may  be  factored  into  a  pair  of  equations  first  order  in  time,  with  the  latter 
instance  giving  the  well  known  Dirac  equation.  In  the  case  of  a  Dirac  particle,  we 
identify  the  frequency  of  probability  of  reversal  a  =  moC^/ih’,  with  moC^  the  rest  mass 
energy  [gaveau].  This  can  provide,  for  example,  an  interpretation  for  the  Zitterbewegung 
phenomena  of  Dirac  theory. 

Classically  or  quantum  mechanically,  lattice  gas  dynamics  may  be  thought  of  in  terms  of 
scattering  due  to  local  potentials.  There  is  an  associated  scattering  matrix,  leading  to 
transmission  and  reflection  coefficients.  Building  upon  such  an  approach,  recent  work 
has  used  quantum  random  walk  to  examine  diffusion  in  ID  crystalline  nanostructures 
[godoy].  The  telegraph  equation  results  in  the  continuum  limit  for  an  irreversible, 
second-order  Markov  process. 

A  QLGA  includes  the  sequential  repetition  of  four  main  steps  [yepez,  yepezOl,  yepez06, 
vahala].  First,  initialization  creates  the  quantum-mechanical  initial  state  that  corresponds 
to  the  intial  probability  distribution  for  a  partial  differential  equation  to  be  solved. 
Secondly,  in  the  collision  step,  a  unitary  transformation  is  applied  in  parallel  to  all  the 
local  Hilbert  spaces  in  the  lattice.  Next,  in  the  measurement  step,  the  quantum  states  of 
all  the  nodes  are  read  out.  Lastly,  these  results  are  shifted  or  “streamed”  to  neighboring 
lattices  sites,  providing  reinitialization  of  the  lattice  in  the  state  which  corresponds  to  the 
updated  probability  distribution. 

QLGAs  have  been  shown  to  solve  the  diffusion.  Burgers,  Boltzmann,  Schroedinger,  and 
Dirac  equations  [yepez,  vahala,  yepezOl,  boghosian98,  yepez02].  In  some  QLGAs  (e.g., 
[boghosian98,  yepez02])  the  measurement  step  is  omitted  and  the  generally  entangled 
quantum  states  are  streamed.  This  places  much  stronger  requirements  on  the  quantum 
computing  hardware,  but  gives  an  exponential  speed  up  over  classical  simulation.  We 
discuss  an  algorithm  with  intermediate-time  measurements. 

There  is  much  applications  interest  in  both  classical  and  QLGAs.  Elsewhere 
[coffeycolburn],  we  have  demonstrated  that  hybrid  computing,  running  versions  of 
diffusion  processing,  could  be  useful  for  the  enhancement  of  digital  images.  In  particular, 
diffusion  of  intensities,  with  a  constraint  on  the  difference  of  pixel  iterate  values,  gives 
selective  smoothing  within  an  image.  In  another  very  recent  scenario,  we  have  developed 
QLGAs  for  the  Maxwell  equations  by  starting  from  a  Dirac  formulation  [coffey].  The 
algorithms  may  be  executed  with  measurement  only  at  the  final  time,  resulting  in  an 
exponential  speed  up  over  classical  simulation. 
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Very  recently,  a  telegraph-diffusion  operator  has  been  proposed  for  purposes  of  image 
restoration  and  denoising  [ratner].  This  approach  requires  the  solution  of  a  nonlinear 
telegraph  equation  with  diffusivity  dependent  upon  the  gradient  of  the  intensity  function. 
Our  results  now  indicate  the  possibility  to  apply  a  QLGA  with  a  signal  strength  constraint 
for  a  ID  telegraph  equation  for  obtaining  the  denoising  of  digital  signals. 

A  drawback  of  the  ordinary  diffusion  Equation  8  for  many  applications  is  that  the 
associated  propagation  speed  is  infinite.  For  any  positive  time  t,  there  can  be  diffusion, 
albeit  usually  very  small,  to  arbitrarily  large  distances.  The  telegraph  equation  offers  one 
way  of  correcting  this  aspect:  it  models  diffusion  with  a  finite  propagation  speed.  This 
can  be  very  important  for  modeling  diffusion  in  a  variety  of  contexts  including  turbulent 
fluids  and  biological  processes  and  ecological  problems  (e.g.,  [codling]).  The  search  for 
a  fully  special  relativistic  diffusion  equation  remains  an  open  and  important  problem  for 
statistical  physics  and  other  areas,  but  the  telegraph  equation  provides  an  improvement 
over  Equation  8. 

Whereas  the  parabolic  Equation  8  has  a  number  of  well  known  properties,  including 
satisfying  a  maximum  principle,  the  behavior  of  solutions  of  the  telegraph  equation  is 
generally  more  complicated.  As  shown  for  instance  by  a  Fourier  series  solution  of  the 
telegraph  equation  with  special  zero  Dirichlet  or  Neumann  boundary  conditions,  there  is 
a  variety  of  behavior  of  the  solutions  of  the  telegraph  equation. 

The  QEGA  has  a  significant  numerical  advantage  inherent  in  its  formulation.  This  is  the 
guaranteed  stability  due  to  the  use  of  a  unitary  collision  operator.  For  a  hyperbolic 
equation  as  we  are  dealing  with  here,  this  is  no  small  matter.  We  recall  that  in 
comparison  an  explicit  finite  difference  scheme  for  a  wave  equation  must  satisfy  the 
Courant-Friedrichs-Eevy  condition  [gasdynbk]  as  a  necessary  constraint.  Roughly 
described,  the  Courant-Friedrichs-Eevy  condition  arises  from  ensuring  that  the  domain  of 
dependence  of  the  numerical  method  contains  the  domain  of  dependence  of  the  partial 
differential  equation  being  solved.  It  has  the  direct  consequence  of  limiting  how  large  the 
time  step  may  be  taken  in  relation  to  the  size  of  the  spatial  discretization.  The  severity  of 
the  Courant-Friedrichs-Eevy  condition  can  be  reduced  only  at  substantial  computational 
cost.  Either  the  time  step  is  drastically  reduced,  or  another  method  such  as  an  implicit 
scheme  is  required.  In  the  latter  event,  there  is  significant  additional  computational  cost 
in  solving  a  set  of  coupled  equations  at  each  time  iteration.  Even  then,  if  the  boundary 
conditions  are  not  treated  fully  implicitly  also,  the  Courant-Friedrichs-Eevy  constraint 
will  become  manifest. 

A  QEGA  also  offers  a  significant  advantage  as  far  as  realizing  a  hybrid  architecture.  This 
is  because  if  the  nodal  qubits  have  sufficiently  long  coherence  time,  no  quantum  error 
correction  is  required.  In  contrast,  many  other  methods  require  quantum  error  correction, 
and  this  is  typically  a  tremendous  increase  in  resource.  Typical  error  correcting 
techniques  encode  one  logical  qubit  in  either  5  or  7  physical  qubits.  On  top  of  this, 
several  levels  of  concatenation  are  used. 


29 


The  ID  telegraph  equation  arises  as  the  probability  density  funetion  (pdf)  for  the 
displaeement  at  time  t  for  persistent  random  walk  on  a  ID  lattice  in  the  continuum  limit. 
With  a  form  of  momentum  introduced  into  the  random  walk,  persistent  random  walk  has 
applications  in  describing  scattering  and  diffusion  in  disordered  media.  As  we  have 
mentioned,  the  telegraph  equation  has  solutions  with  a  finite  velocity  of  propagation,  and 
this  can  provide  an  advantage  for  describing  heat  propagation,  light  dispersion  in  turbid 
media,  or  in  biological  modeling. 

We  have  discussed  a  new  QLGA  for  the  ID  telegraph  equation  that  subsumes  one  for  the 
diffusion  equation,  while  complementing  that  for  the  ID  Dirac  equation.  Both  the 
telegraph  and  Dirac  equations  may  be  developed  from  microscopic  models  with  particles 
undergoing  random  reversals  of  direction.  The  resulting  QLGA  for  the  ID  telegraph 
equation  is  highly  parallelizable  and  well  founded  on  physical  principles  including 
conservation  laws.  This  algorithm  offers  the  prospect  for  simulation  on  combined 
classical-quantum  computing  architectures.  In  such  a  computing  environment,  local 
nodes  with  two  qubits  each  are  connected  to  nearest  neighbors  with  classical 
communication. 

We  have  verified  our  QLGA  on  model  problems  which  take  into  account  the  boundary 
conditions  and  the  initial  conditions  on  both  the  solution  and  its  first  order  time  derivative 
[cc2009].  Test  problems  included  those  from  references  [goldstein]  and  [morette]. 
Despite  using  an  approximation  to  implement  non  periodic  boundary  conditions  for  some 
test  problems,  we  were  still  able  to  find  accurate  solutions.  The  incorporation  of  further 
types  of  boundary  conditions  is  an  area  for  future  research. 
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4.  Summary 


The  concepts  of  entropy  and  of  diffusion  have  been  common  themes  throughout  much  of 
this  research.  Diffusion  and  the  accompanying  idea  of  scale  space  find  many  applications 
in  image  sharpening  and  restoration  and  other  processing  tasks.  Entropy  measures 
provide  a  useful  monitor  of  diffusion  processing,  as  well  as  a  way  to  understand  the 
information  content  of  images. 

In  the  overall  scope  of  the  investigations,  we  were  able  to  make  contributions  to  quantum 
information  science  over  a  fairly  broad  range  of  topics,  from  a  new  high  level  algorithm 
to  detailed  logic  gate  decompositions.  The  vast  majority  of  these  topics  are  additionally 
complemented  with  external  publications. 

As  a  byproduct,  we  made  various  contributions  to  the  analytic  evaluation  of  reduced 
Feynman  diagram  integrals.  In  doing  so  we  advanced  the  theory  of  some  particular 
special  functions  of  mathematical  physics  and  elaborated  the  theory  of  generalized 
harmonic  number  sums. 

We  found  often  that  much  ‘ground  work’  needs  to  be  done — there  is  an  absence  of 
quantum  algorithmic  ‘infrastructure’.  Although  quantum  computing  is  a  very  active  area, 
just  which  applications  are  most  amenable  to  quantum  approaches  often  remain  elusive. 
Two  example  questions  from  the  latter  portion  of  our  overall  investigation  include; 
Which  scheduling  and  logistics  problems  have  oracle  functions  that  may  be  implemented 
efficiently  in  order  to  perform  adaptive  search?  Which  type  of  limited  logistics  problem 
may  yet  have  a  polynomial-time  quantum  solution? 

Lastly  we  made  specific  contributions  to  quantum  lattice  gas  algorithms.  We  showed 
how  the  important  Maxwell  equations  of  electromagnetism  could  be  solved  with  this 
class  of  algorithms.  We  further  generalized  an  algorithm  for  the  diffusion  equation  so 
that  the  more  complicated  telegraph  equation  could  be  simulated  with  a  QLGA.  The 
latter  partial  differential  equation  has  multiple  applications,  including  transport  in 
disordered  media,  image  enhancement,  heat  propagation,  and  elsewhere. 

As  a  result  of  this  project,  senior  student  Ben  James  Jones  received  the  Ryan  Sayers 
memorial  award  for  research  as  a  combined  computer  science  and  engineering  physics 
major. 
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